#ifndef UTILITY_H
#define UTILITY_H
#include <stdio.h>
#include <stdlib.h>
#include <list>
#include <cctype>
using namespace std;

extern "C" void hello(void) {
	printf("Hello, world!\n");
}

//****************************************************************************
template<typename T>
void PrintVectorFixedChars(std::ostream& ostr, const std::vector<T>& myvec,
		const std::string& delimiter, int length) {
	//****************************************************************************

	vector < string > vec;

	for (int k = 0; k < myvec.size(); k++) {
		string shortned;
		shortned = (myvec[k]);
		shortned.resize(length, ' ');
		vec.push_back(shortned);
	}

	std::copy(vec.begin(), vec.end(),
			std::ostream_iterator < T > (ostr, delimiter.c_str()));
}
// use:	PrintVector(std::cout, v, " ");

//****************************************************************************
template<typename T>
T round(T input, float unit) {
	//****************************************************************************
	T res;

	float remnant;
	int division;
	division = floor(input / unit);

	remnant = input - division * unit;

	if (remnant / unit <= 0.5) {

		res = division;
	} else if (remnant / unit > 0.5) {
		res = division + unit;
	}

	return res;

}
// use:	PrintVector(std::cout, v, " ");

//****************************************************************************
template<typename T>
void PrintVectorFixedDigits(std::ostream& ostr, const std::vector<T>& vec,
		const std::string& delimiter) {
	//****************************************************************************

	ostr << fixed;

	std::copy(vec.begin(), vec.end(),
			std::ostream_iterator < T > (ostr, delimiter.c_str()));
	ostr.unsetf(ios_base::floatfield);

}
// use:	PrintVector(std::cout, v, " ");

//****************************************************************************
template<typename T>
void PrintVector(std::ostream& ostr, const std::vector<T>& vec,
		const std::string& delimiter) {
	//****************************************************************************
	std::copy(vec.begin(), vec.end(),
			std::ostream_iterator < T > (ostr, delimiter.c_str()));
}
// use:	PrintVector(std::cout, v, " ");

//****************************************************************************
void PrintColumns(std::ostream& efficiencies_table_file,
		vector<vector<string> > efficiencies_table_column) {
	//****************************************************************************

	for (int riga = 0; riga < efficiencies_table_column[0].size(); riga++) {
		for (int colonna = 0; colonna < efficiencies_table_column.size();
				colonna++) {
			//efficiencies_table_file << fixed;
			efficiencies_table_file.unsetf(ios::floatfield);
			//efficiencies_table_file << setprecision(3);
			efficiencies_table_file.precision(3);
			efficiencies_table_file << setprecision(3)
					<< efficiencies_table_column[colonna][riga] << '\t';
		}
		efficiencies_table_file << endl;
		efficiencies_table_file.unsetf(ios_base::floatfield);

	}
}

//****************************************************************************
bool repeting_particle(vector<int> vector1, vector<int> vector2) {
	//****************************************************************************
	for (int pos1 = 0; pos1 < vector1.size(); pos1++) {
		for (int pos2 = 0; pos2 < vector2.size(); pos2++) {
			if (vector1[pos1] == vector2[pos2]) {
				return true;
			}
		}
	}
	return false;
}

//****************************************************************************
float charge(int pid) {
	//****************************************************************************
	float ch = 9999;
	if (pid == 21 || pid == 22 || pid == 23 || abs(pid) == 12 || abs(pid) == 14
			|| abs(pid) == 16) {

		ch = 0.;
	}

	if (pid == 11 || pid == 13 || pid == 15 || pid == -23) {
		ch = -1.;
	}
	if (pid == -11 || pid == -13 || pid == -15 || pid == 23) {
		ch = 1.;
	}

	if (pid == 1 || pid == 3 || pid == 5) {
		ch = -1. / 3;
	}

	if (pid == -1 || pid == -3 || pid == -5) {
		ch = 1. / 3;
	}

	if (pid == 2 || pid == 4 || pid == 6) {
		ch = 2. / 3;
		//cout << pid << "->" << ch << endl;

	}

	if (pid == -2 || pid == -4 || pid == -6) {
		ch = -2. / 3;
	}

	return ch;

}

//****************************************************************************
void one_digit(float f) {
	//****************************************************************************
	cout << setprecision(1) << f;
	cout << setprecision(5);
}

//****************************************************************************
void printvector(vector<int> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<double> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<bool> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<pair<int, int> > v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i].first << " " << v[i].second;
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<long double> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<float> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printvector(vector<string> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << " ";
		i++;
	}
	cout << endl;
}
template<typename T>
//****************************************************************************
void printvectornl(vector<T> v) {
	//****************************************************************************
	int i = 0;
	while (i < v.size()) {
		cout << v[i] << endl;
		i++;
	}

}

//****************************************************************************
float avgmatrix(vector<vector<float> > matrix) {
	//****************************************************************************
	float result = 0;
	int riga = 0;
	int colonna = 0;
	int elements = 0;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[0].size()) {
			elements++;
			result = result + matrix[riga][colonna];
			colonna++;
		}
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return result / elements;
}

//****************************************************************************
TLorentzVector zBoost(TLorentzVector FV, long double Y) {
	//****************************************************************************
	TLorentzVector BoostedFV;
	long double Yprime;
	Yprime = FV.Rapidity() + Y;
	BoostedFV.SetE(FV.Et() * cosh(Yprime));
	BoostedFV.SetPz(FV.Et() * sinh(Yprime));
	BoostedFV.SetPx(FV.Px());
	BoostedFV.SetPy(FV.Py());
	return BoostedFV;
}

//****************************************************************************
TLorentzVector BoostZ(TLorentzVector FV, long double Y) {
	//****************************************************************************
	TLorentzVector BoostedFV;
	//	long double Yprime;
	//	Yprime = FV.Rapidity() + Y;
	BoostedFV.SetE(FV.E() * cosh(Y) - FV.Pz() * sinh(Y));
	BoostedFV.SetPz(-FV.E() * sinh(Y) + FV.Pz() * cosh(Y));
	BoostedFV.SetPx(FV.Px());
	BoostedFV.SetPy(FV.Py());
	return BoostedFV;
}

//****************************************************************************
float avgcolumn(vector<vector<float> > matrix, int column) {
	//****************************************************************************
	float result = 0;
	int riga = 0;
	int elements = 0;
	while (riga < matrix.size()) {
		elements++;
		result = result + matrix[riga][column];
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return result / elements;
}

//****************************************************************************
float variancematrix(vector<vector<float> > matrix) {
	//****************************************************************************
	float avg;
	avg = avgmatrix(matrix);
	float result = 0;
	int riga = 0;
	int colonna = 0;
	int elements = 0;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[0].size()) {
			elements++;
			result = result + pow((matrix[riga][colonna] - avg), 2);
			colonna++;
		}
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return sqrt(result / elements);
}

//****************************************************************************
float varianceofacolumn(vector<vector<float> > matrix, int column) {
	//****************************************************************************
	float avg;
	avg = avgcolumn(matrix, column);
	float result = 0;
	int riga = 0;
	//int colonna=0;
	int elements = 0;
	while (riga < matrix.size()) {
		elements++;
		result = result + pow((matrix[riga][column] - avg), 2);
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return sqrt(result / elements);
}

//****************************************************************************
int elementsofamatrix(vector<vector<float> > matrix, float lowerbound,
		float upperbound) {
	//****************************************************************************

	float result = 0;
	int riga = 0;
	int colonna = 0;
	int elements = 0;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[0].size()) {

			if (matrix[riga][colonna] < upperbound
					&& matrix[riga][colonna] > lowerbound) {
				elements++;
			}

			colonna++;
		}
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return elements;
}

//****************************************************************************
int elementsofamatrix(vector<vector<float> > matrix) {
	//****************************************************************************

	int riga = 0;
	int colonna = 0;
	int elements = 0;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[0].size()) {
			elements++;
		}
		riga++;
	}
	//cout<<result<<endl;
	//cout<<"elems "<<elements<<endl;
	return elements;
}

//****************************************************************************
int maxvector(vector<int> matrix) {
	//****************************************************************************
	int result = matrix[0];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = max(result, matrix[riga]);

		riga++;
	}
	return result;
}

//****************************************************************************
float maxvector(vector<float> matrix) {
	//****************************************************************************
	float result = matrix[0];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = max(result, matrix[riga]);

		riga++;
	}
	return result;
}

//****************************************************************************
float minvector(vector<float> matrix) {
	//****************************************************************************
	float result = matrix[0];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = min(result, matrix[riga]);

		riga++;
	}
	return result;
}
//****************************************************************************
int minvector(vector<int> matrix) {
	//****************************************************************************
	int result = matrix[0];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = min(result, matrix[riga]);

		riga++;
	}
	return result;
}

//****************************************************************************
vector<float> constant_times_vector(float c, vector<bool> A) {
	//****************************************************************************
	vector<float> result;

	int i = 0;

	//cout<<"riga "<<i<<endl;

	while (i < A.size()) {
		//cout<<"colonna "<<j<<endl;
		result.push_back(c * 1.0 * A[i]);
		i++;
	}

	return result;
}

//****************************************************************************
vector<float> constant_times_vector(float c, vector<float> A) {
	//****************************************************************************
	vector<float> result;

	int i = 0;

	//cout<<"riga "<<i<<endl;

	while (i < A.size()) {
		//cout<<"colonna "<<j<<endl;
		result.push_back(c * 1.0 * A[i]);
		i++;
	}

	return result;
}

//****************************************************************************
vector<float> constant_times_vector(float c, vector<double> A) {
	//****************************************************************************
	vector<float> result;

	int i = 0;

	//cout<<"riga "<<i<<endl;

	while (i < A.size()) {
		//cout<<"colonna "<<j<<endl;
		result.push_back(c * 1.0 * A[i]);
		i++;
	}

	return result;
}

//****************************************************************************
float maxcolumnofamatrix(vector<vector<float> > matrix, int column) {
	//****************************************************************************
	float result = matrix[0][column];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = max(result, matrix[riga][column]);

		riga++;
	}
	return result;
}

//****************************************************************************
float mincolumnofamatrix(vector<vector<float> > matrix, int column) {
	//****************************************************************************
	float result = matrix[0][column];
	int riga = 0;

	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {

		result = min(result, matrix[riga][column]);

		riga++;
	}
	return result;
}

//****************************************************************************
float maxmatrix(vector<vector<float> > matrix) {
	//****************************************************************************
	float result = matrix[0][0];
	int riga = 0;
	int colonna = 0;
	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[riga].size()) {

			if (matrix[riga][colonna] > result) {
				//cout << "at position " << riga << "," << colonna;
				//cout << " was " << result;
				result = matrix[riga][colonna];
				//cout << " now is " << result << endl;
			}

			//result = max(result, matrix[riga][colonna]);

			colonna++;
		}
		riga++;
	}

	return result;
}

//****************************************************************************
float minmatrix(vector<vector<float> > matrix) {
	//****************************************************************************
	float result = matrix[0][0];
	int riga = 0;
	int colonna = 0;
	//cout<<matrix.size()<<endl;
	//cout<<matrix[0].size()<<endl;
	while (riga < matrix.size()) {
		colonna = 0;
		while (colonna < matrix[riga].size()) {
			//cout<<riga<<colonna<<endl;
			//cout<< matrix[riga][colonna]<<endl;
			result = min(result, matrix[riga][colonna]);
			colonna++;
		}
		riga++;
	}
	return result;
}

//****************************************************************************
vector<float> progressive_efficiencies_to_steps(vector<float> v) {
	//****************************************************************************

	vector<float> step;

	step.push_back(1.);

	//cout<<"righe "<<v.size()<<endl;
//	for (int i = v.size() - 1; i >= 0; i--) {
	for (int i = 0; i < v.size() - 1; i++) {

		step.push_back(v[i + 1] / v[i]);

	}

	return step;
}

//****************************************************************************
string vector2CSV(vector<int> v) {
	//****************************************************************************

	string strout;

	//cout<<"righe "<<v.size()<<endl;
	for (int i = 0; i < v.size(); i++) {
		char thisitem[5];
		sprintf(thisitem, "%i", v[i]);

		strout.append(thisitem);
		strout.append(",");

	}
	string::iterator itend;
	itend = strout.end() - 1;
	strout.erase(itend);

	return strout;
}

//****************************************************************************
string matrix2string(vector<vector<int> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;

	string strout;

	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			char thisitem[5];
			sprintf(thisitem, "%i", v[i][j]);

			strout.append(thisitem);
			strout.append(",");

			j++;
		}
		string::iterator itend;
		itend = strout.end() - 1;
		strout.erase(itend);
		strout.append(";");
		i++;
	}

	return strout;
}

template<typename T>
//****************************************************************************
void printmatrixgen(vector<vector<T> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			cout << " " << v[i][j] << " ";
			j++;
		}
		cout << endl;
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printmatrix(vector<vector<string> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			cout << " " << v[i][j] << " ";
			j++;
		}
		cout << endl;
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printmatrix(vector<vector<bool> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			cout << " " << v[i][j] << " ";
			j++;
		}
		cout << endl;
		i++;
	}
	cout << endl;
}

//****************************************************************************
void printmatrix(vector<vector<int> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			cout << " " << v[i][j] << " ";
			j++;
		}
		cout << endl;
		i++;
	}
	cout << endl;
}

//****************************************************************************
vector<pair<int, int> > matrixelements_less_than_one(vector<vector<float> > v) {
	//****************************************************************************
	vector < pair<int, int> > result;

	int i = 0;
	int j = 0;

	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			if (v[i][j] < 1) {
				pair<int, int> elem;
				elem.first = i;
				elem.second = j;
				result.push_back(elem);
			}
			j++;
		}

		i++;
	}
	return result;
}

//****************************************************************************
vector<vector<float> > initialize_matrix(float value, int rows, int columns) {
	//****************************************************************************
	vector < vector<float> > A;
	if (A.empty() == true) {
		int i = 0;
		int j = 0;
		//cout<<"righe "<<v.size()<<endl;
		while (i < rows) {
			//cout<<"riga "<<i<<endl;
			//cout<<"di dimensione"<<v[i].size()<<endl;
			vector<float> firstlevel_result;
			j = 0;
			while (j < columns) {
				//cout<<"colonna "<<j<<endl;
				firstlevel_result.push_back(value);
				j++;
			}

			A.push_back(firstlevel_result);
			i++;
		}
	} else if (A.empty() == false) {
		cout << "the matrix has already been initialized" << endl;
	}
	return A;
}

//****************************************************************************
vector<float> initialize_vector(float value, int rows) {
	//****************************************************************************
	vector<float> A;
	if (A.empty() == true) {
		int i = 0;
		//int j = 0;
		//cout<<"righe "<<v.size()<<endl;
		while (i < rows) {
			//cout<<"riga "<<i<<endl;
			//cout<<"di dimensione"<<v[i].size()<<endl;
			A.push_back(value);
			i++;
		}
	} else if (A.empty() == false) {
		cout << "the matrix has already been initialized" << endl;
	}
	return A;
}

////****************************************************************************
//vector<float> initialize_vector_array(const int input[]) {
//	//****************************************************************************
//
//	//int myints[] = input;
//
//	vector<int> A(input, input + sizeof(input) / sizeof(int));
//
//	return A;
//}

//****************************************************************************
int element_minimum_vector(vector<long double> distance) {
	//****************************************************************************

	long double mindistance;
	int best_k, res;
	mindistance = distance[0];
	best_k = 0;

	for (int k = 0; k < distance.size(); k++) {
		//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

		if (1. / distance[k] > 1. / mindistance) {
			best_k = k;
			mindistance = distance[k];
		}

	} //for

	res = best_k;

	return res;

} //best k find

//****************************************************************************
void printmatrix(vector<vector<float> > v) {
	//****************************************************************************
	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < v.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		j = 0;
		while (j < v[i].size()) {
			//cout<<"colonna "<<j<<endl;
			cout << " " << v[i][j] << " ";
			j++;
		}
		cout << endl;
		i++;
	}
	cout << endl;
}

//****************************************************************************
vector<vector<float> > summatrix(vector<vector<float> > B,
		vector<vector<float> > A) {

	//****************************************************************************
	vector < vector<float> > result;

	if (A.size() != B.size()) {
		cout << "sum of matrices of different size is not possible" << endl;
		exit(1);
	}
	for (int elem = 0; elem < B.size(); elem++) {
		if (A[elem].size() != B[elem].size()) {
			cout << "sum of matrices of different size is not possible" << endl;
			exit(1);
		}
	}

	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		vector<float> firstlevel_result;
		j = 0;
		while (j < A[i].size()) {
			//cout<<"colonna "<<j<<endl;
			firstlevel_result.push_back(A[i][j] + B[i][j]);
			j++;
		}

		result.push_back(firstlevel_result);
		i++;
	}

	return result;
}

//****************************************************************************
template<typename T>
int sign(T number) {

	//****************************************************************************
	if (number > 0) {
		return 1;
	}

	if (number < 0) {
		return -1;
	}

}

//****************************************************************************
template<typename T>
int any2closestinteger(T number) {

	//****************************************************************************
	T c, f, m;
	f = floor(number);
	c = ceil(number);

	m = min(c - number, number - f);
	if (m == c - number) {
		int cint;
		cint = (int) c;
		return cint;
	}
	if (m == number - f) {
		int fint;
		fint = (int) f;
		return fint;
	}

}
//****************************************************************************
template<typename T>
vector<T> mergevector(vector<T> A, vector<T> B) {

	//****************************************************************************
	vector<T> result;

	int i = 0;

	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		//cout<<"colonna "<<j<<endl;
		result.push_back(A[i]);
		i++;
	}
	i = 0;
	while (i < B.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		//cout<<"colonna "<<j<<endl;
		result.push_back(B[i]);
		i++;
	}

	return result;
}

//****************************************************************************
vector<float> sumvector(vector<float> B, vector<float> A) {

	//****************************************************************************
	vector<float> result;

	if (A.size() != B.size()) {
		cout << "sum of matrices of different size is not possible" << endl;
		exit(1);
	}

	int i = 0;

	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		//cout<<"colonna "<<j<<endl;
		result.push_back(A[i] + B[i]);
		i++;
	}

	return result;
}

//****************************************************************************
vector<vector<float> > constant_times_matrix(float c,
		vector<vector<float> > A) {

	//****************************************************************************
	vector < vector<float> > result;

	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		vector<float> firstlevel_result;
		j = 0;
		while (j < A[i].size()) {
			//cout<<"colonna "<<j<<endl;
			firstlevel_result.push_back(c * A[i][j]);
			j++;
		}

		result.push_back(firstlevel_result);
		i++;
	}

	return result;
}

//****************************************************************************
vector<vector<float> > matrix1overmatrix2(vector<vector<float> > A,
		vector<vector<float> > B) {

	//****************************************************************************
	vector < vector<float> > result;

	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		vector<float> firstlevel_result;
		j = 0;
		while (j < A[i].size()) {
			//cout<<"colonna "<<j<<endl;
			firstlevel_result.push_back(A[i][j] / B[i][j]);
			j++;
		}

		result.push_back(firstlevel_result);
		i++;
	}

	return result;
}

//****************************************************************************
vector<vector<float> > sqrt_of_matrix(vector<vector<float> > A) {

	//****************************************************************************
	vector < vector<float> > result;

	int i = 0;
	int j = 0;
	//cout<<"righe "<<v.size()<<endl;
	while (i < A.size()) {
		//cout<<"riga "<<i<<endl;
		//cout<<"di dimensione"<<v[i].size()<<endl;
		vector<float> firstlevel_result;
		j = 0;
		while (j < A[i].size()) {
			//cout<<"colonna "<<j<<endl;
			firstlevel_result.push_back(sqrt(A[i][j]));
			j++;
		}

		result.push_back(firstlevel_result);
		i++;
	}

	return result;
}

//****************************************************************************
//****************** GENERATE THE SUBSETS ************************************
//****************************************************************************
void ksub_next(int n, int k, vector<int> & a, bool *more)
//****************************************************************************
//****************************************************************************
//****************************************************************************

//
//  Purpose:
//
//    KSUB_NEXT generates the subsets of size K from a set of size N.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    29 May 2003
//
//  Author:
//
//    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms for Computers and Calculators,
//    Second Edition,
//    Academic Press, 1978,
//    ISBN: 0-12-519260-6,
//    LC: QA164.N54.
//
//  Parameters:
//
//    Input, int N, the size of the set from which subsets are drawn.
//
//    Input, int K, the desired size of the subsets.  K must
//    be between 0 and N.
//
//    Output, int A[K].  A[I] is the I-th element of the
//    subset.  Thus A[I] will be an integer between 1 and N.
//    Note that the routine will return the values in A
//    in sorted order: 1 <= A[0] < A[1] < ... < A[K-1] <= N
//
//    Input/output, bool *MORE.  Set MORE = FALSE before first call
//    for a new sequence of subsets.  It then is set and remains
//    TRUE as long as the subset computed on this call is not the
//    final one.  When the final subset is computed, MORE is set to
//    FALSE as a signal that the computation is done.
//
		{

	int j;
	static int m = 0;
	static int m2 = 0;

	if (k < 0 || n < k) {
		cout << "\n";
		cout << "KSUB_NEXT - Fatal error!\n";
		cout << "N = " << n << "\n";
		cout << "K = " << k << "\n";
		cout << "but 0 <= K <= N is required!\n";
		exit(1);
	}

	if (!(*more)) {
		m2 = 0;
		m = k;

//		hello();
//	    cout << k << "out of "<< n<< endl;
	} else {
		if (m2 < n - m) {
			m = 0;
		}
		m = m + 1;
		m2 = a[k - m];
	}

	for (j = 1; j <= m; j++) {
		//cout<<k<<j<<m<<endl;
		//cout<<a.empty();
		a.insert(a.begin() + k + j - m - 1, m2 + j);

	}

	*more = (a[0] != (n - k + 1));

//	cout << "more in ksubnext is "<< *more << endl;

	return;
}

//****************************************************************************
void moveparticle(int mostenergeticparticle, vector<int> & notused,
		vector<int> & jet) {
	//****************************************************************************

	vector<int>::iterator it;

	jet.push_back(mostenergeticparticle);

	it = find(notused.begin(), notused.end(), mostenergeticparticle);
	notused.erase(it);

}

/* generic functio to flatten an any-dim vector to a 1 dim vector*/
/*
 //****************************************************************************
 template<typename T>
 vector<T> flatten2vector(const std::vector<T>& vec) {
 //****************************************************************************
 T object;
 object = vec[0];

 for (int level = 1; object.size() != 1; level++) {
 object = object[0];
 }

 //now object is the innermost quantity and we can extract its type

 typedef typeinner
 vector<T> output;

 for (int elem = 0; elem < vec.size(); elem++) {

 if (elem != elemtoskip) {

 output.push_back(vec[elem]);
 }

 }

 return output;

 }*/

//****************************************************************************
template<typename T>
vector<T> matrix2vector(const std::vector<vector<T> > & vec) {
	//****************************************************************************
	//std::copy(vec.begin(), vec.end(), std::ostream_iterator<T>(ostr, delimiter.c_str()));

	vector<T> output;

	for (int riga = 0; riga < vec.size(); riga++) {

		for (int colonna = 0; colonna < vec[riga].size(); colonna++) {

			output.push_back(vec[riga][colonna]);
		}

	}

	return output;

}

template<typename T>
//****************************************************************************
vector<T> stripthroughlist(const std::vector<T> & kplets,
		vector<bool> kplet_pass_thecut) {
	//****************************************************************************

	vector<T> result;
	result.resize(0);

	if (kplets.size() == kplet_pass_thecut.size()) {
		/* sono almeno di stessa dimensione! */
		for (int kpletcounter = 0; kpletcounter < kplet_pass_thecut.size();
				kpletcounter++) {

			if (kplet_pass_thecut[kpletcounter] == true) {
				result.push_back(kplets[kpletcounter]);
			}

		}

	} else {
		puts("the list of successes and the list of kplets do not match");
		exit(1);
	}

	return result;

}
template<typename T>
//****************************************************************************
bool compare_pairs_VF(pair<long double, T> p1, pair<long double, T> p2) {
	//****************************************************************************
	long double pT1 = p1.first;
	long double pT2 = p2.first;

	if (pT1 < pT2) {
		return true;
	} else {
		return false;
	}
}

//****************************************************************************
vector<vector<int> > homemade_pair_sort(
		vector<pair<long double, vector<int> > > inputpairs) {
	//****************************************************************************
	vector<vector<int> >::iterator it;

	vector<long double> values;
	vector<long double>::iterator values_it;

	if (inputpairs.size() == 1) {

		vector < vector<int> > result;

		values.push_back(inputpairs[0].first);
		result.push_back(inputpairs[0].second);

		return result;

	} else if (inputpairs.size() == 2) {

		vector < vector<int> > result;

		if (inputpairs[0].first >= inputpairs[1].first) {
			values.push_back(inputpairs[0].first);
			result.push_back(inputpairs[0].second);
			values.push_back(inputpairs[1].first);
			result.push_back(inputpairs[1].second);
		} else if (inputpairs[0].first < inputpairs[1].first) {
			values.push_back(inputpairs[1].first);
			result.push_back(inputpairs[1].second);
			values.push_back(inputpairs[0].first);
			result.push_back(inputpairs[0].second);
		}

		return result;

	} else if (inputpairs.size() > 2) {

		vector < vector<int> > result;

		if (inputpairs[0].first >= inputpairs[1].first) {
			values.push_back(inputpairs[0].first);
			result.push_back(inputpairs[0].second);
			values.push_back(inputpairs[1].first);
			result.push_back(inputpairs[1].second);
		} else if (inputpairs[0].first < inputpairs[1].first) {
			values.push_back(inputpairs[1].first);
			result.push_back(inputpairs[1].second);
			values.push_back(inputpairs[0].first);
			result.push_back(inputpairs[0].second);
		}

		//now fill before the first encoutered value that is less thant the value at hand
		//find the position where to fill (the iterator!)

		//cout << "dopo due elementi il vettore e`" << endl;
		//printvector(values);
		//cout << endl;
		//	cout<<"dopo due elementi il vettore e`"<<endl;printmatrix(result);cout<<endl;

		for (int kvalue = 2; kvalue < inputpairs.size(); kvalue++) {
			it = result.begin();
			for (values_it = values.begin();
					values_it < values.end()
							&& inputpairs[kvalue].first < *values_it;
					values_it++) {
				//cout<<"valori allo step di"<< ""<<":"<<inputpairs[kvalue].first <<" "<<*values_it<<endl;
				it++;
			}
			//		cout << "the search stopped at the elment " << *it;
			// this will stop when the value at hand, inputpairs[kvalue].first, is larger that the object in the list with iterator it
			// I will insert the value just before this element
			vector<int> tempvint;
			tempvint = inputpairs[kvalue].second;

			result.insert(it, tempvint);
			//add element to values to match with result

			long double tempvalue;
			tempvalue = inputpairs[kvalue].first;
			values.insert(values_it, tempvalue);
			//cout << "dopo " << kvalue << "elementi il vettore dei kplets e`" << endl;
			//printmatrix(result);
			//cout << endl;

			//cout << "dopo " << kvalue << "elementi il vettore dei valori e`" << endl;
			//printvector(values);
			//cout << endl;

		}

		return result;
	}

}

//template<typename T>
////****************************************************************************
//vector<T> homemade_pair_sorttemp(vector<pair<long double, T> > inputpairs) {
//	//****************************************************************************
//	typedef typename std::vector<T>::const_iterator cIt;
//	std::vector<vector<vector<int> > >::const_iterator it;
//
//	vector<long double> values;
//	vector<long double>::iterator values_it;
//
//	if (inputpairs.size() == 1) {
//		cout<<"one value to sort"<<endl;
//		vector<T> result;
//		values.push_back(inputpairs[0].first);
//		result.push_back(inputpairs[0].second);
//		return result;
//	}
//
//	if (inputpairs.size() == 2) {
//		cout<<"two values to sort"<<endl;
//		vector<T> result;
//		if (inputpairs[0].first > inputpairs[1].first) {
//			values.push_back(inputpairs[0].first);
//			result.push_back(inputpairs[0].second);
//			values.push_back(inputpairs[1].first);
//			result.push_back(inputpairs[1].second);
//		} else if (inputpairs[0].first < inputpairs[1].first) {
//			values.push_back(inputpairs[1].first);
//			result.push_back(inputpairs[1].second);
//			values.push_back(inputpairs[0].first);
//			result.push_back(inputpairs[0].second);
//		}
//
//		return result;
//	}
//
//	if (inputpairs.size() > 2) {
//		cout<<"more than 2 values to sort"<<endl;
//		vector<T> result;
//		if (inputpairs[0].first > inputpairs[1].first) {
//			values.push_back(inputpairs[0].first);
//			result.push_back(inputpairs[0].second);
//			values.push_back(inputpairs[1].first);
//			result.push_back(inputpairs[1].second);
//		} else if (inputpairs[0].first < inputpairs[1].first) {
//			values.push_back(inputpairs[1].first);
//			result.push_back(inputpairs[1].second);
//			values.push_back(inputpairs[0].first);
//			result.push_back(inputpairs[0].second);
//		}
//
//		//now fill before the first encoutered value that is less thant the value at hand
//		//find the position where to fill (the iterator!)
//		cout<<"dopo due elementi il vettore e`"<<endl;printvector(result);cout<<endl;
//
//		for (int kvalue = 2; kvalue < inputpairs.size(); kvalue++) {
//			cIt it = result.begin();
//			for (values_it = values.begin(); values_it < values.end()
//					&& inputpairs[kvalue].first < *values_it; values_it++) {
//				it++;
//			}
//			// this will stop when the value at hand, inputpairs[kvalue].first, is larger that the object in the list with iterator it
//			// I will insert the value just before this element
//			//			T tempvalue;
//			//			tempvalue=inputpairs[kvalue].second;
//			//			result.insert(it,tempvalue);
//
//		}
//
//		return result;
//	}
//
//}

template<typename T>
//****************************************************************************
vector<T> sort2ndaccordingto1st(vector<long double> object_to_sort,
		vector<T> previousclassmemeber) {
	//****************************************************************************
	//make pairs and then sort them

	vector < pair<long double, T> > pairs_to_sort;

	//	typedef typename std::list<pair<long double,T> >::const_iterator cIt;

	if (object_to_sort.size() == previousclassmemeber.size()) {

		for (int obj = 0; obj < object_to_sort.size(); obj++) {
			pair<long double, T> temppair;
			temppair.first = object_to_sort[obj];
			temppair.second = previousclassmemeber[obj];
			pairs_to_sort.push_back(temppair);
		}

		vector<T> result;

		//		pairs_to_sort.sort(compare_pairs_VF);
		//		for (cIt it=pairs_to_sort.begin(); it!=pairs_to_sort.end(); ++it){
		//
		//			result.push_back(it->second);
		//		}

		result = homemade_pair_sort(pairs_to_sort);

		return result;

	}
}

template<typename T>
//****************************************************************************
vector<T> sort2ndaccordingto1st(vector<double> object_to_sort,
		vector<T> previousclassmemeber) {
	//****************************************************************************
	//make pairs and then sort them

	vector < pair<long double, T> > pairs_to_sort;

	//	typedef typename std::list<pair<long double,T> >::const_iterator cIt;

	if (object_to_sort.size() == previousclassmemeber.size()) {

		for (int obj = 0; obj < object_to_sort.size(); obj++) {
			pair<long double, T> temppair;
			temppair.first = object_to_sort[obj];
			temppair.second = previousclassmemeber[obj];
			pairs_to_sort.push_back(temppair);
		}

		vector<T> result;

		//		pairs_to_sort.sort(compare_pairs_VF);
		//		for (cIt it=pairs_to_sort.begin(); it!=pairs_to_sort.end(); ++it){
		//
		//			result.push_back(it->second);
		//		}

		result = homemade_pair_sort(pairs_to_sort);

		return result;

	}
}

template<typename T>
//****************************************************************************
vector<T> sort2ndaccordingto1st(vector<float> object_to_sort,
		vector<T> previousclassmemeber) {
	//****************************************************************************
	//make pairs and then sort them

	vector < pair<long double, T> > pairs_to_sort;

	//	typedef typename std::list<pair<long double,T> >::const_iterator cIt;

	if (object_to_sort.size() == previousclassmemeber.size()) {

		for (int obj = 0; obj < object_to_sort.size(); obj++) {
			pair<long double, T> temppair;
			temppair.first = object_to_sort[obj];
			temppair.second = previousclassmemeber[obj];
			pairs_to_sort.push_back(temppair);
		}

		vector<T> result;

		//debug
//		for (int elem = 0; elem < pairs_to_sort.size(); elem++) {
//
//			long double primo = pairs_to_sort[elem].first;
//			T secondo = pairs_to_sort[elem].second;
//
//			cout << "(" << primo << ", ";
//			PrintVector(std::cout, secondo, " ");
//			cout << " )" << endl;
//
//		}

		result = homemade_pair_sort(pairs_to_sort);

		return result;

	}
}

//****************************************************************************
template<typename T>
vector<T> stripelement(int elemtoskip, const std::vector<T>& vec) {
	//****************************************************************************
	//std::copy(vec.begin(), vec.end(), std::ostream_iterator<T>(ostr, delimiter.c_str()));

	vector<T> output;

	for (int elem = 0; elem < vec.size(); elem++) {

		if (elem != elemtoskip) {

			output.push_back(vec[elem]);
		}

	}

	return output;

}

/********************************************/
string rootparticle2lhe(int partid, vector<TRootLHEFParticle*> partvector)
/********************************************/
{
	char output[] = "";

	char PID[] = "";
	char Status[] = "";
	char Mother1[] = "";
	char Mother2[] = "";
	char ColorLine1[] = "";
	char ColorLine2[] = "";
	char Px[] = "";
	char Py[] = "";
	char Pz[] = "";
	char Mass[] = "";
	char E[] = "";

	sprintf(PID, "%d ", partvector[partid]->PID);
	sprintf(Status, "%d ", partvector[partid]->Status);
	sprintf(Mother1, "%d ", partvector[partid]->Mother1);
	sprintf(Mother2, "%d ", partvector[partid]->Mother2);
	sprintf(ColorLine1, "%d ", partvector[partid]->ColorLine1);
	sprintf(ColorLine2, "%d ", partvector[partid]->ColorLine2);
	sprintf(Px, "%f ", partvector[partid]->Px);
	sprintf(Py, "%f ", partvector[partid]->Py);
	sprintf(Pz, "%f ", partvector[partid]->Pz);
	sprintf(E, "%f ", partvector[partid]->E);
	sprintf(Mass, "%f ", partvector[partid]->M);

	strcat(output, PID);
	strcat(output, Status);
	strcat(output, Mother1);
	strcat(output, Mother2);
	strcat(output, ColorLine1);
	strcat(output, ColorLine2);
	strcat(output, Px);

	strcat(output, Py);

	strcat(output, Pz);

	strcat(output, E);
	strcat(output, Mass);
	strcat(output, "0 9");

	return output;

}

//****************************************************************************
void findkplets(vector<vector<vector<int> > > & allkplets, int k,
		int nmembers) {
	//****************************************************************************

	vector < vector<int> > kplets;
	//alla fine multiplet va dentro allkplets

	bool *more = new bool;
	*more = false;
	vector<int> out;
	//find the first multiplet with ksub_next(SET,SIZE_OF_SUBSET, OUTPUT, FLAG_TO_CONTINUE)

	//cout << "asking for " << k << " out of " <<nmembers << " in findkplets"<< endl;
	//cout << "the bool is "<< *more << endl;

	ksub_next(nmembers, k, out, more);

	//cout << "the bool is "<< *more << endl;

	kplets.push_back(out);
//	cout<<"so far the results is"<<endl;
//	printvector(out);

	//find the remaining multiplets with ksub_next()
	while (*more) {
		ksub_next(nmembers, k, out, more);
		out.resize(k);
		kplets.push_back(out);
		//cout<<"so far the results is"<<endl;
		//printvector(out);
	}
	//cout<<"the kplets are"<<endl;
	//printmatrix(kplets);
	//fill the container of the multiplets
	allkplets.push_back(kplets);
} //EndOf findkplets(allkplets,k,n)

//****************************************************************************
void largeventfindkplets(vector<vector<vector<int> > > & allkplets, int k,
		int nmembers) {
	//****************************************************************************

	vector < vector<int> > kplets;
	//alla fine multiplet va dentro allkplets

	bool *more = new bool;
	vector<int> out;

	if (k < 3 || nmembers - k < 3) {

		//cout<<"k was "<<k<<" and we will use it!"<<endl;

		//find the first multiplet with ksub_next()
		ksub_next(nmembers, k, out, more);
		kplets.push_back(out);

		//find the remaining multiplets with ksub_next()
		while (*more) {
			ksub_next(nmembers, k, out, more);
			out.resize(k);
			kplets.push_back(out);
		}

		//printmatrix(kplets);
		//fill the container of the multiplets
		allkplets.push_back(kplets);

	} else {
		//cout<<"k was "<<k<<" and we will NOT use it!"<<endl;
		vector < vector<int> > merenda;
		vector<int> merendina;
		merendina.push_back(100000000);
		merenda.push_back(merendina);

		allkplets.push_back(merenda);
	}

} //EndOf findkplets(allkplets,k,n)

//****************************************************************************
bool isinrange(int x, int low, int up) {
	//****************************************************************************

	if (x <= up && x >= low) {
		return 1;
	} else if (x > up || x < low) {
		return 0;
	}
}

//****************************************************************************
bool isinrange(float x, float low, float up) {
	//****************************************************************************

	if (x <= up && x >= low) {
		return 1;
	} else if (x > up || x < low) {
		return 0;
	}
}

typedef long double (*FUN)(vector<vector<int> >, vector<TLorentzVector>,
		vector<TRootLHEFParticle*>);

//****************************************************************************
bool meets_the_requirement(float results, string minmaxrange, float OBSmin,
		float OBSmax) {
	//****************************************************************************
	if (minmaxrange == "min") {
		if (results > OBSmin) {
			return true;
		} else {
			return false;
		}

	} else if (minmaxrange == "max") {

		if (results < OBSmax) {
			return true;
		} else {
			return false;
		}

	} else if (minmaxrange == "range") {

		if ((results > OBSmin) && (results < OBSmax)) {
			return true;
		} else {
			return false;
		}

	}

}

//*************************************************
bool vectorcontains(vector<bool> cuts_results, bool query) {
	//*************************************************

	for (int i = 0; i < cuts_results.size(); i++) {
		if (cuts_results[i] == query) {
			return 1;
		}
	}
	return 0;
}

template<typename T>
//*************************************************
bool vectorcontains(vector<T> list, T query) {
	//*************************************************

	for (int i = 0; i < list.size(); i++) {
		if (list[i] == query) {
			return true;
		}
	}
	return false;
}

template<typename T>
//*************************************************
pair<bool, int> vectorcontainsatposition(vector<T> list, T query) {
	//*************************************************

	for (int i = 0; i < list.size(); i++) {
		if (list[i] == query) {

			pair<bool, int> output;
			output.first = true;
			output.second = i;
			return output;
		}
	}

	pair<bool, int> output;
	output.first = false;
	output.second = -1;
	return output;
}

//*************************************************
vector<int> get_CSV(char inputchars[]) {
	//*************************************************
	vector<int> result;
	istringstream instream;
	string inputstring;
	inputstring = inputchars;
	instream.str(inputstring);
	string tempstring;
	int tempint;

	while (getline(instream, tempstring, ',')) {
		result.push_back(atoi(tempstring.c_str()));
	}
	return result;

}

//*************************************************
string int2string(int f) {
	//*************************************************
	ostringstream buffer;
	buffer << f;
	string str = buffer.str();
	return str;
}

//*************************************************
vector<string> vector2vectorstring(vector<float> f) {
	//*************************************************
	vector < string > res;
	for (int i = 0; i < f.size(); i++) {
		ostringstream buffer;
		buffer << fixed;
		buffer << f[i];
		res.push_back(buffer.str());
		buffer.unsetf(ios_base::floatfield);

	}
	return res;
}

//*************************************************
vector<string> vector2vectorstringprecision(vector<float> f) {
	//*************************************************
	vector < string > res;
	for (int i = 0; i < f.size(); i++) {
		ostringstream buffer;
		buffer.precision(3);
		buffer << f[i];
		res.push_back(buffer.str());
		buffer.unsetf(ios_base::floatfield);

	}
	return res;
}

//*************************************************
string vector2string(vector<int> f) {
	//*************************************************
	ostringstream buffer;
	for (int i = 0; i < f.size(); i++) {
		buffer << f[i] << ",";
	}
	string str = buffer.str();
	return str;
}

//*************************************************
string vector2stringprecision(vector<int> f) {
	//*************************************************
	ostringstream buffer;
	buffer.precision(3);
	for (int i = 0; i < f.size(); i++) {
		buffer << f[i] << ",";
	}
	string str = buffer.str();
	return str;
}

//*************************************************
vector<string> uniform_length_string(vector<string> vec) {
	//*************************************************
	vector < string > res;
	int max_length = 0;
	for (int i = 0; i < vec.size(); i++) {

		if (vec[i].size() > max_length) {
			max_length = vec[i].size();
		}
	}

	res = vec;

	for (int i = 0; i < res.size(); i++) {

		res[i].resize(max_length, ' ');

	}

	return res;
}

//*************************************************
string float2string(float f) {
	//*************************************************
	ostringstream buffer;
	buffer << f;
	string str = buffer.str();
	return str;
}

//*************************************************
vector<int> get_separated_values(string inputstring, char separator) {
	//*************************************************
	vector<int> result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;
	int tempint;
	//	char[1000] charseparator="";
	//
	//	strcat(charseparator, separator.c_str() );
	//
	//	char charseparator = (char)separator.c_str();
	//	char sep=separator&;
	while (getline(instream, tempstring, separator)) {
		result.push_back(atoi(tempstring.c_str()));
	}
	return result;

}

bool is_integer(float k) {
	return std::floor(k) == k;
}

//*************************************************
vector<string> get_separated_strings(string inputstring, char separator) {
	//*************************************************
	vector < string > result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;
	int tempint;
	//	char[1000] charseparator="";
	//
	//	strcat(charseparator, separator.c_str() );
	//
	//	char charseparator = (char)separator.c_str();
	//	char sep=separator&;
	while (getline(instream, tempstring, separator)) {
		result.push_back(tempstring);
	}
	return result;

}

//*************************************************
vector<int> get_CSV(string inputstring) {
	//*************************************************
	vector<int> result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;
	int tempint;

	while (getline(instream, tempstring, ',')) {
		result.push_back(atoi(tempstring.c_str()));
	}
	return result;

}

//*************************************************
map<string, float> assignment2map(vector<string> assignments) {
	//*************************************************
	map<string, float> mymap;

	for (int kass = 0; kass < assignments.size(); kass++) {
		vector < string > fields;
		string tempfield;

		//		cout << kass << " " << assignments[kass] << endl;

		istringstream instream;
		instream.str(assignments[kass]);

		while (getline(instream, tempfield, '=')) {
			//cout<<"tempfield is " <<tempfield<<endl;
			if (tempfield != "" && tempfield != " ") {
				fields.push_back(tempfield);
				//cout<<"filling fields with "<<tempfield<<endl;
			}
		}
		mymap.insert(pair<string, float>(fields[0], atof(fields[1].c_str())));
		//		cout << fields[0] << " " << fields[1] << endl;
	}

	return mymap;
}

//*************************************************
map<string, string> assignment2mapS(vector<string> assignments) {
	//*************************************************

	map < string, string > mymap;

	for (int kass = 0; kass < assignments.size(); kass++) {
		vector < string > fields;
		string tempfield;
		istringstream instream;
		instream.str(assignments[kass]);

		while (getline(instream, tempfield, '=')) {
			//cout<<"tempfield is " <<tempfield<<endl;
			if (tempfield != "" && tempfield != " ") {
				fields.push_back(tempfield);
				//cout<<"filling fields with "<<tempfield<<endl;
			}
		}
		mymap.insert(pair<string, string>(fields[0], fields[1]));
	}

	return mymap;
}

//*************************************************
vector<float> get_CSVf(string inputstring) {
	//*************************************************
	vector<float> result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;
	int tempint;

	while (getline(instream, tempstring, ',')) {
		result.push_back(atof(tempstring.c_str()));
	}
	return result;

}
////*************************************************
//float * vector2array(vector<float> s) {
//	//*************************************************
//	int n;
//	n = s.size();
//	float result[n];
//
//	for (int k = 0; k < s.size(); k++) {
//		result[k] = s[k];
//	}
//
//	return result;
//}

//*************************************************
float get_ValueFromEquation(string inputstring) {
	//*************************************************
	float result;
	vector < string > fields;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;
	int tempint;

	while (getline(instream, tempstring, '=')) {
		fields.push_back(tempstring);
	}
	if (fields.size() == 2) {
		result = atof(fields[1].c_str());
		return result;
	} else {
		cout << inputstring
				<< " is not a valid string to interpret for valorization "
				<< endl;
		exit(1);
	}

}

//*********************************************
bool is_a_parton(int PID) {
	//*********************************************

	if ((abs(PID) <= 6) || (PID == 21)) {
		return true;
	} else {
		return false;
	}

}

//****************************************************************************
vector<vector<int> > oneplets2eventnumbering(vector<vector<int> > allkplets,
		vector<vector<int> > classmembers) {
	//****************************************************************************
	vector < vector<int> > result;
	int colonna = 0;

	for (int riga = 0; riga < allkplets.size(); riga++) {

		result.push_back(classmembers[allkplets[riga][colonna] - 1]);
		// the content of the first level of class members is already a vector
	}

	return result;
}

//*************************************************
bool isinstring(string s1, string s2) {
	//*************************************************
	bool result = false;
	if (s1.find(s2) != std::string::npos) {
		result = true;
	}
	return result;
}

//*************************************************
bool isinclass(int PidN, vector<int> classpids) {
	//*************************************************
	int i = 0;
	while (i < classpids.size()) {
		if (PidN == classpids[i]) {
			return true;
		}
		if (classpids[i] == 0) {
			cout << " all PIDs are taken" << endl;
			return true;
		}
		i++;
	}
	return false;
}

template<typename T>
//*************************************************
void debugprint(string str, T val) {
	//*************************************************
	cout << str << val << endl;
}

//*************************************************
bool isinclass(vector<int> classpids, int PidN) {
	//*************************************************
	int i = 0;
	while (i < classpids.size()) {
		if (PidN == classpids[i]) {
			return true;
		}
		if (classpids[i] == 0) {
			cout << " all PIDs are taken" << endl;
			return true;
		}
		i++;
	}
	return false;
}

//*************************************************
bool isinvector(vector<int> v, int k) {
	//*************************************************
	int i = 0;
	while (i < v.size()) {
		if (k == v[i]) {
			return true;
		}

		i++;
	}
	return false;
}

//*************************************************
pair<string, string> get_function_name_and_arguments(string inputstring) {
	//*************************************************
	pair < string, string > name_and_args;

	istringstream instream, instreamargs;
	instream.str(inputstring);
	string tempstring;
	vector < string > fields;
	while (getline(instream, tempstring, '[')) {
		//cout << tempstring << endl;
		fields.push_back(tempstring);
	}
	string args;
	instreamargs.str(fields[1]);
	getline(instreamargs, args, ']');

	name_and_args.first = fields[0];
	name_and_args.second = args;

	return name_and_args;
}

//*************************************************
vector<string> getCSstrings(string inputstring) {
	//*************************************************
	vector < string > result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;

	while (getline(instream, tempstring, ',')) {
		//cout << tempstring << endl;
		result.push_back(tempstring);
	}
	return result;
}
//*************************************************
vector<string> get_USstrings(string inputstring) {
	//*************************************************
	vector < string > result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;

	while (getline(instream, tempstring, 'U')) {
		//cout << tempstring << endl;
		result.push_back(tempstring);
	}
	return result;
}
//*************************************************
vector<string> get_pipeSstrings(string inputstring) {
	//*************************************************
	vector < string > result;
	istringstream instream;
	instream.str(inputstring);
	string tempstring;

	while (getline(instream, tempstring, '|')) {
		//cout << tempstring << endl;
		result.push_back(tempstring);
	}
	return result;
}

//****************************************************************************
int howmany(vector<bool> v, bool value) {
	//****************************************************************************
	int counter = 0;
	for (int i = 0; i < v.size(); i++) {
		if (v[i] == value) {
			counter++;
		}
	}
	return counter;
}

//****************************************************************************
int howmanyint(vector<int> v, int value) {
	//****************************************************************************
	int counter = 0;
	for (int i = 0; i < v.size(); i++) {
		if (v[i] == value) {
			counter++;
		}
	}
	return counter;
}

//****************************************************************************
void valorize(int & tobevalorized, string name,
		map<string, float> global_parameters_map) {
	//****************************************************************************
	map<string, float>::iterator it;
	it = global_parameters_map.find(name);
	if (it != global_parameters_map.end()) {
		tobevalorized = global_parameters_map.find(name)->second;
	} else if (it == global_parameters_map.end()) {
		cout << "warning: " << name
				<< " has no assigned value and will take the default (if any)"
				<< endl;

	}
}

//****************************************************************************
void valorize(float & tobevalorized, string name,
		map<string, float> global_parameters_map) {
	//****************************************************************************
	map<string, float>::iterator it;
	it = global_parameters_map.find(name);
	if (it != global_parameters_map.end()) {
		tobevalorized = global_parameters_map.find(name)->second;
	} else if (it == global_parameters_map.end()) {
		cout << "warning: " << name
				<< " has no assigned value and will take the default (if any)"
				<< endl;

	}
}

//****************************************************************************
void valorizeS(string & tobevalorized, string name,
		map<string, string> global_parameters_map) {
	//****************************************************************************
	map<string, string>::iterator it;
	it = global_parameters_map.find(name);
	if (it != global_parameters_map.end()) {
		tobevalorized = global_parameters_map.find(name)->second;
	} else if (it == global_parameters_map.end()) {
		cout << "warning: " << name
				<< " has no assigned value and will take the default (if any)"
				<< endl;

	}
}

//****************************************************************************
bool findeventnum(int eventnum, Int_t & rootindex, TFile *f, string treename,
		string branchname) {
	//****************************************************************************

	//	f = TFile::Open(jetfilename, "READ");

	//	TTree *t1 = (TTree*) f->Get("tvec");

	char treenamechar[100] = "";
	strcat(treenamechar, treename.c_str());

	char branchnamechar[100] = "";
	strcat(branchnamechar, branchname.c_str());

	TTree *t1 = (TTree*) f->Get(treenamechar);

	Int_t ev;

	TBranch *bvpx = 0;
	t1->SetBranchAddress("eventnumber", &ev, &bvpx);

	//	cout << "cerco l'evento #" << eventnum << endl;

	//	t1->SetBranchAddress("eventnumber", &ev);

	//read all entries and fill the histograms
	Int_t nentries = (Int_t) t1->GetEntries();

	//	cout << "ho trovato un numero di entries " << nentries << endl;

	for (Int_t i = 0; i < nentries; i++) {
		//		hello();

		Long64_t tentry = t1->LoadTree(i);
		bvpx->GetEntry(tentry);

		//		t1->GetEntry(i);

		if (ev == eventnum) {

			rootindex = i;
			//			f->Close("R");
			//			delete f;
			t1->ResetBranchAddresses();
			return true;
		}

	}
	t1->ResetBranchAddresses();
	return false;

}

//****************************************************************************
void printlv(TLorentzVector v) {
	//****************************************************************************
	//	std::cout << v[0] << " " << v[1] << " " << v[2] << " " << v[3] << std::endl;
	std::cout << v.E() << " " << v.Px() << " " << v.Py() << " " << v.Pz()
			<< std::endl;

}
//****************************************************************************
void print3v(TVector3 bCoM) {
	//****************************************************************************

	cout << bCoM.X() << " " << bCoM.Y() << " " << bCoM.Z() << endl;
}
//****************************************************************************
TVector3 vectorproduct(TLorentzVector A, TLorentzVector B) {
	//****************************************************************************
	TVector3 result;
	long double paz, pay, pax, pbx, pby, pbz;

	pax = A.Px();
	pay = A.Py();
	paz = A.Pz();
	pbx = B.Px();
	pby = B.Py();
	pbz = B.Pz();

	result.SetX(pay * pbz - paz * pby);
	result.SetY(paz * pbx - pax * pbz);
	result.SetZ(pax * pby - pay * pbx);

	return result;

}

//****************************************************************************
vector<vector<int> > getthepairs(int ntobeclulstered) {
	//****************************************************************************

	//		cout << "the pairing into " << ntobeclulstered << " elements has to be generated" << endl;
	vector < vector<vector<int> > > allkplets;
	findkplets(allkplets, 2, ntobeclulstered);

	return allkplets[0];

}

//****************************************************************************
vector<vector<int> > getthepairsfromfile(int ntobeclulstered, TFile *filepair) {
	//****************************************************************************

	vector < vector<int> > returnpairings;

	TTree *treepairs;
	filepair->GetObject("pairs", treepairs);

	bool pairing_exists;
	pairing_exists = false;

	Int_t rootindex;

	//	hello();
	pairing_exists = findeventnum(ntobeclulstered, rootindex, filepair, "pairs",
			"pairings");

	if (pairing_exists == false) {

		//		cout << "the pairing into " << ntobeclulstered << " elements has to be generated" << endl;
		vector < vector<vector<int> > > allkplets;
		findkplets(allkplets, 2, ntobeclulstered);
		returnpairings = allkplets[0];

		return returnpairings;

	} else if (pairing_exists == true) {

		//		cout << "the pairing for N="<< ntobeclulstered<<" will be reused" << endl;

		std::vector < std::vector<int> > *pairings = 0;
		TBranch *bpairs = 0;
		treepairs->SetBranchAddress("pairings", &pairings, &bpairs);

		Long64_t tentry = treepairs->LoadTree(rootindex);
		bpairs->GetEntry(tentry);

		returnpairings = *pairings;

		//		filepair->Close("R");
		//		delete filepair;

		return returnpairings;

		//				printmatrix(allpairs);

	}
}

//****************************************************************************
bool CheckOneDistinct(vector<vector<int> > kpletsthatpassthecut,
		int cuthowmany) {
	//****************************************************************************

	/* faccio i cuthowmanypletti */
	vector < vector<vector<int> > > chplets;
	chplets.resize(0);

	int nmember;
	nmember = kpletsthatpassthecut.size();
	//cout<<cuthowmany+1<< " "<<nmember<<endl;
	// I add plus 1 because the request > and not >= is assumed
	findkplets(chplets, cuthowmany + 1, nmember);

	/* costruisco le liste flat con i c.h.plets */

	/* loop sulle liste*/
	for (int listcounter = 0; listcounter < chplets[0].size(); listcounter++) {
		vector<int> currentlist;
		currentlist.resize(0);

		for (int i = 0; i < chplets[0][listcounter].size(); i++) {

			currentlist = mergevector(currentlist,
					kpletsthatpassthecut[chplets[0][listcounter][i] - 1]);

		}

		/* controllo che nessun elemento si ripeta*/

		// cout<<"checking if this is distinct"<<endl;
		//printvector(currentlist);
		bool accept_this_list;
		accept_this_list = true;
		for (int listelement = 0; listelement < currentlist.size();
				listelement++) {
			if (howmanyint(currentlist, currentlist[listelement]) > 1) {
				accept_this_list = false;
				// cout<<"is NOT distinct"<<endl;

			}
		}

		if (accept_this_list == true) {
			// cout<<"is distinct"<<endl;

			return true;
		}

	}

	return false;

}

//****************************************************************************
int pid23timescharge(int pid) {
	//****************************************************************************
	if ((pid == 1) || (pid == 3) || (pid == 5)) {

		return -1;
	}

	if ((pid == -1) || (pid == -3) || (pid == -5)) {
		return 1;
	}

	if ((pid == -2) || (pid == -4) || (pid == -6)) {
		return -2;
	}

	if ((pid == 2) || (pid == 4) || (pid == 6)) {

		//		hello();
		return 2;
	}

	if (pid == 11 || pid == 13 || pid == 15 || pid == 24) {
		return 3;
	}

	if (pid == -11 || pid == -13 || pid == -15 || pid == -24) {
		return -3;
	}
	if (pid == 12 || pid == 14 || pid == 16 || pid == -12 || pid == -14
			|| pid == -16 || pid == 21 || pid == 22 || pid == 23) {
		return 0;
	} else {
		cout << "unknown charge for PID " << pid << endl;
		return 1000;

	}

}

//*************************************************
float get_float_parameter(string parameters, string searched_parameter) {
	//*************************************************
	/* FIND THE FUNCTION FOR THE OBSERVABLE*/
	vector < string > tempassignments = getCSstrings(parameters);
	map<string, float> tempmymap;
	tempmymap = assignment2map(tempassignments);
	// from here on the variables are initialized and valorized from the string provided as input
	float Fparameter;
	valorize(Fparameter, searched_parameter, tempmymap);

	return Fparameter;

}

//*************************************************
string get_string_parameter(string parameters, string searched_parameter) {
	//*************************************************
	/* FIND THE FUNCTION FOR THE OBSERVABLE*/
	vector < string > tempassignments = getCSstrings(parameters);
	map < string, string > tempmymap;
	tempmymap = assignment2mapS(tempassignments);
	// from here on the variables are initialized and valorized from the string provided as input
	string Sparameter;
	valorizeS(Sparameter, searched_parameter, tempmymap);

	return Sparameter;

}

//*************************************************
FUN get_the_obs_function(string parameters) {
	//*************************************************
	/* FIND THE FUNCTION FOR THE OBSERVABLE*/
	vector < string > tempassignments = getCSstrings(parameters);
	map < string, string > tempmymap;
	tempmymap = assignment2mapS(tempassignments);
	// from here on the variables are initialized and valorized from the string provided as input
	string OBS;
	valorizeS(OBS, "OBS", tempmymap);
	//	string minmaxrange;
	//	valorizeS(minmaxrange, "minmaxrange", tempmymap);
	void *dl_handle;
	//Dynamically open myself
	dl_handle = dlopen(NULL, RTLD_LAZY);

	if (!dl_handle) {
		cout << "Failed to open myself" << endl;
		exit(1);
	}
	char *thisobs_funcname;
	thisobs_funcname = const_cast<char *>(OBS.c_str());
	//Look up the function the user asked for
FUN 	thisobs_func = (FUN) dlsym(dl_handle, thisobs_funcname);
	if (!thisobs_func) {
		cout << " Failed to find the desired obervable " << OBS
				<< "for the generic filter " << endl;
		exit(1);
	}

	return thisobs_func;

}
/* END OF FIND THE FUNCTION FOR THE OBSERVABLE*/

//*************************************************
int MaxElementInAVector(vector<float> pTsJ) {
	//*************************************************

	int maxelem;
	float max;
	max = pTsJ[0];
	maxelem = 0;

	for (int i = 0; i < pTsJ.size(); i++) {

		if (pTsJ[i] > max) {

			max = pTsJ[i];
			maxelem = i;

		}

	}

	return maxelem;
}

//*************************************************
int MaxElementInAVectorWithExceptions(vector<float> pTsJ,
		vector<int> exceptions) {
	//*************************************************

	int maxelem;
	float max;

	for (int i = 0; i < pTsJ.size(); i++) {

		if (!(isinvector(exceptions, i))) {

			max = pTsJ[i];
			maxelem = i;

		}

	}

	for (int i = 0; i < pTsJ.size(); i++) {

		if (pTsJ[i] > max && !(isinvector(exceptions, i))) {

			max = pTsJ[i];
			maxelem = i;

		}

	}

	return maxelem;
}

//*************************************************
vector<int> MaxKElementsInAVector(vector<float> pTsJ, int k) {
//************************************************

	vector<int> res;

	res.resize(0);

	vector<float> input;

	input = pTsJ;

	if (k <= input.size()) {

		for (int i = 0; i < k; i++) {
			//	cout << "The input is ";
			//	PrintVector(std::cout, input, " ");
			//	cout << endl;

			int thismaxelem;
			thismaxelem = MaxElementInAVectorWithExceptions(input, res);

//			cout << " the max element is " << thismaxelem << endl;

			res.push_back(thismaxelem);

			//		cout << "The result is ";
			//	PrintVector(std::cout, res, " ");
			//cout << endl;

			//	input.erase(input.begin() + thismaxelem);

		}

	} else {
		cout << "Too few elements: cannot find the " << k
				<< " largest elements in  a set of " << pTsJ.size() << endl;
	}

	return res;

}
//*************************************************

#endif
